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Abstract 

Food scientists at the U.S. Army's Natick Solider Center have devel- 
oped a model for the lifecyle of the bacteria Staphylococcus aureus in 
| intermediate moisture bread. In this article, we study this model us- 

ing dynamical systems and Lie symmetry methods. We calculate center 
manifolds and Lie symmetries for different cases of parameter values and 
compare our results to those of the food scientists. 



1 Introduction 



1.1 The model 

To ensure the food that U.S. soldiers receive is as safe as possible, the growth 
of bacteria such as Staphylococcus aureus (S. aureus) needs to be addressed. 
The system of equations considered in this paper arises from a "quasi-chemical" 
kinetics model for the phases of the microbial life cycle of S. aureus in intermedi- 
ate moisture bread. The food scientists who developed the model confirmed its 
usefulness by fitting to it the data from observations on bread crumbs with vary- 
ing conditions of water activity, pH and temperature. The model differs from 
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previous models in its attempt to model continuously the growth and death of 
the microorganism rather than focusing solely on either growth or inactivation. 
The model was developed by food scientists at the Natick Soldier Center, see 
Taub, et al [4] for more on their techniques. 

The model arose from the observation of four phases in the life cycle of S. 
aureus. The cells pass through the various stages of metabolizing (M), mul- 
tiplying (M*), sensitization to death (M**), and dead (D). Additionally, the 
scientists hypothesized that there was an antagonist (A) present that would af- 
fect the cells. They found that without this added element their original model 
did not fit the observed data with any accuracy. The first step in the process 
describes cells moving from lag phase to growth phase (M — » M*). In the next 
step, cells multiply via binary division and then the newly multiplied cells in- 
teract with an antagonist (M* — » 2M* + A). The last two steps represent two 
different pathways to death: the first with cells interacting with an antagonist, 
then passing to sensitization before death (A + M* — > M** — > D) and lastly the 
cells experiencing natural death (M* — » D). 

The following equations represent the velocities of each of the above steps 
(v) as they relate to the concentrations of cells in various the phases. Each 
equation has a rate constant (k) associated to it. 



vi = kiM (1) 

v 2 = k 2 M* (2) 

v 3 = (W- 9 )k 3 M*A (3) 

Vi = kiM* (4) 

Finally these velocities are represented by the following system of ordinary dif- 
ferential equations: 

M = -vi = -kiM (5) 

M* = vi + v 2 -v 3 -vi = kiM + M* (G - eA) (6) 

A = v 2 - v 3 = M*(k 2 -eA) (7) 

D = v 3 + vi = M*(k 4 + eA) (8) 



where G = k 2 — ki is the net natural growth rate and e = lCP 9 /^. It is assumed 
that all the rate constants have non-negative values. The initial conditions at 
time zero are M(0) = /, the inoculum level I w 10 3 — 10 4 , and M*(0) = A(0) — 
D(0) = 0. 

1.2 A simplification 

We notice that the fourth equation is uncoupled since there are no terms involv- 
ing the variable D in any of the other equations and D depends on M* and A. 
Therefore to investigate the dynamics of our system, we reduce to a system of 
three equations. Renaming our variables (yi = M, 1/2 = M*, y 3 = A) we have 
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(9) 



2 Normalizing the system 

To consider the invariant manifold structure of a system, it is necessary to write 
the system in normal form as follows: 

x = Ax + g(x,y) (10) 
V = By + j(x,y) (11) 

with (x,y) S R™ x R m , the n x n matrix A having eigenvalues with zero real 
part and the m x m matrix B having eigenvalues with nonzero real part. The 
functions g(x,y) and j(x,y) must be zero with zero first partial derivatives at 
the origin. 

The system (|9]) above is not in normal form since the y' 3 equation corresponds 
to the zero eigenvalue piece and the nonlinear term of y' 3 does not have all 
zero partial derivatives at the origin. Thus we must normalize by a change 
of coordinates using the eigenvectors of the matrix of the linear terms of the 
equation. We will investigate the invariant manifolds in a neighborhood of 
G = 0. Writing our system in normal form for nonzero G does not depend on 
the sign of G, so we treat the negative and positive case simultaneously. We 
let T be the matrix of eigenvectors of the eigenvalues of the linear terms of our 
system and let 



(12) 



Using the inverse of the matrix T we can solve for u, v and w, find their 
derivatives and finally write our system in normal form as follows: 

= 0-u+(l-^]f(u,v,w) (13) 

(14) 
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where f(u, v, w) = —-^(v — kiw)(k2V + G(u + k2w)). Since f(u, v, w) and its first 
partials with respect to u, v, and w are all zero at the origin (it, v, w) — (0, 0, 0), 
we have our system in normal form and we see immediately that we have a 
one-dimensional center manifold in the case that G ^ 0. For G > 0, we also 
have a one-dimensional stable and a one-dimensional unstable manifold. For 
G < 0, we have a two-dimensional stable manifold. The system reduced to the 
center manifold simply becomes 



u' = 0. 



(15) 
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For the case G = 0, we have a slightly simpler system of equations: 

(16) 
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Note that we now have two zero eigenvalues and one negative eigenvalue for 
the matrix in the linear term. Since zero is a repeated eigenvalue, we must 
use generalized eigenvectors to find the normalization of this system. Three 
such eigenvectors are (0,0, 1), (0, 1,0), and (1, —1, |^). Then to transform our 
system we again let T be the matrix consisting of these eigenvectors and let 



(17) 



As above, this allows us to write our system in normal form: 

k 2 \ f u\ f g(u,v,w) \ , . 

J {v ) + { g(u,v,w) J [l *> 

w' = -k lW + (19) 

where g(u, v, w) = —e(v — w)(u+ j^w). Since g(u, v, w) is zero at the origin and 
all of its first partial derivatives are also zero at the origin, we can see that we 
have a two dimensional center manifold and a one dimensional stable manifold. 



3 Center manifold calculations 

Recall that a center manifold W c — {(x, y)\y — h{x)} is described by h(x) where 
h(0) — Dh(0) — 0. We consider a system written in normal form 

x = Ax + g(x,y) (20) 
y = By + j(x,y) (21) 

with A having eigenvalues with zero real part and B eigenvalues with nonzero 
real part. Then we determine h(x) by finding the function that satisfies the 
following condition: 

(Mh)(x) = Dh(x)[Ax + g(x,h(x))] - Bh(x) -j(x,h{x)) =0. 

The sign of G does not change the outcome of this calculation, thus we treat 
the case G ^ at once. We have h : V — > M 2 , V C K a neighborhood of the 
origin. Thus let h(x) = (h 1 (x),h 2 (x)) = {ax 2 + bx 3 + 0(x 4 ), cx 2 + dx 3 + 0{x A )). 
Then f(x, hi(x), h 2 (x)) = e(— a + kic)x 3 + 0(x 4 ) resulting in 



(Mh)(x) = 
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Solving for (Mh)(x) — 0, h\(x) = h 2 (x) = 0(x 4 ). Thus up to third order, we 
have h\(x) = h 2 (x) — 0, so a center manifold is simply the u— axis. 

Next we consider the case G = 0. Here h : V — > M, V C R 2 , a neighborhood 
of the origin. We let = h(xi,x 2 ) = ax\ + bx% + cx\x 2 + dx^ + ex\ + Jx\x2 + 
jx\x\. Then we calculate 

/, f .x/ \ / 7 / \ u f \ \ ( k 2X2 + g(xi,X2,h(xi,X 2 )) \ . , N 

(M/0(*) = (M*i>*a). h X2 (x u x 2 ))-^ g{xuX2MxuX2)) )+hh(x u x 2 ) 

— (2ak 2 + ck\)x\X2 + {ak\)x\ + (ck2 + bki)x 2 + {dk\)x\ + [jk 2 + ek\)x\ 
+ (3dk 2 - 2ae - ce + fki)xjx 2 + (2/fe 2 — ce — 2be + jki)xix\ 

resulting in h{x\,x 2 ) = 0(x i ), thus h(xi,x 2 ) — up to order three. Hence in 
this case the uv— plane is a center manifold. 

4 Lie Symmetry 

Recall that a Lie symmetry is a map from the set of solutions of a system 
of differential equations to the set itself. For a system of first order ordinary 
differential equations 

y' k =ujk(t,yi,y2,---,yn), k=l,...,n (22) 

the Lie symmetries that transform the variables t, y±, . . . , y n have infinitesimal 
generators of the form 

X = id t + md Vl + md V2 H h Vnd Vn (23) 

where £ = f (t, yi,y 2 ,..., y n ) and r] k = rj k (t, yi,y 2 , ■ ■■,y n ) for all k. The in- 
finitesimal generator must satisfy the Linearized Symmetry Condition: 

X^(y' k -cj k )=0, k = l,...,n (24) 

when (f2"l?l) holds. In this case the prolongation of X is as follows: 

X« = X + + + • • • + V ^d v , n (25) 

where is defined as 77^ = — y' k D t t;. The total derivative D t in this 

case is D t — dt + y[d yi + • • • + y' n d Vn . Thus we have the following: 

% = dtVk+y[d Vl rik+y 2 dy 2 r}k-\ \-y'n d ynVk-y'k(dt^+y'id yi ^+y 2 d.y 2 ^-\ \-y' n d y ^). 

(26) 

A system of first order ODEs has an infinite number of symmetries. We 
find symmetries by solving for the functions £, r\ k that satisfy the Linearized 
Symmetry Condition (|24|) . This condition reduces to a system of PDEs which 
are computationally difficult to solve. We use the "Intro to Symmetry" package 
in Mathematica and a script included in Cantwell [1] to calculate the symmetries 
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for our system. We are limited in the symmetries we can calculate by our 
computing power. In the case we calculate symmetries up to third order 
in our original coordinates t/i , y 2 , and y 3 and then use a change of coordinates 
on our symmetries to rewrite in the coordinates u, v, and w of our equations in 
normal form. Since the case G = involves simpler equations, we are able to 
calculate these symmetries directly from the equations in normal form, however 
we followed the same method as in the G^O case since we want to be able to 
compare cases. 



4.1 The case G^O 

The infinitesimals of the Lie symmetries (up to order 3) are listed in an array 
with {£, f]i, r]2, 773}, representing the infinitesimal generator X = £d t + Vid yi + 

X x = {1,0,0,0} 

X2 = {y2,~k 1 y 1 y 2 ,k 1 y 1 y 2 + Gyl - £yly 3 ,k 2 yl - eyjy 3 } 
X3 = {y3, -hym, kiyiys + Gy 2 y 3 - ey 2 yl, k 2 y 2 y 3 - £2/21/3} 

r „ G e k 2 e 

Xa = {0, -2/1,2/1 + —?/2 - -^-2/22/3, y x m ~ y l V2V ^ 

. 1 G e k 2 e . 

^5 = {-j^t, -ty\,ty\ + —ty 2 - ^2/22/3, -^2/2 - j^ty 2 y 3 } 

— 1 ki 2 ~ki 2 G , ~k 2 , 
a 6 = { — 2/i,—2/i, 2/i 2/12/2 + 2/12/22/3, 2/12/2 + 2/12/22/3} 

£ £ £ £ £ 

Then we transform the infinitesimal generators of the Lie symmetries found in 
the yi coordinates as follows. If X is an infinitesimal generator in y i: then X = 
(Xt)dt + (Xu)d u + (Xv)d v + (Xw)d w is the corresponding infinitesimal generator 
for a Lie symmetry in the u,v,w coordinates The transformed symmetries in 
the form X = {£, rf 1 , rf 2 , 1J3} where X = £d t + 771 <9„ + if 2 d v + rf 3 d w : 

X x = {1,0,0,0} 

X 2 = {j(u, v, w), q{G~ k 2 )j{u, v, w)f(u, v, w),j(u, v, w)(Gv + f(u, v, w)), -kiwj(u, v, w)} 
X 3 = {l(u, v, w), t;(G — k 2 )l(u, v, w)f(u, v, w),l(u, v, w)(Gv + f(u, v, to)), —k\wl{u, v, to)} 

Gr 

Xt = {0, -^—(G-k 2 )f(u, v, to), ^-(Gv + f(u,v, w)),-w } 
Gfci Ki 

X 5 = {^,-t-(G - k 2 )f(u,v,w),-^(Gv + f(u,v,w)),-tw} 
fci Gfci fci 

X 6 = {m(u, v, to), — k2)m(u, v, w)f(u, v, w),m(u, v, w)(Gv + f(u, v, to)), —k\wm(u, v, to) 

G 

where f(u, v, to) is as above, j(u, v,w) = v — kiw, l(u, w, to) = u + ^jV + k2tu 
and m(u, v,w) — ~-(G + ki)w. 



CENTER MANIFOLDS AND LIE SYMMETRIES 



7 



4.2 The case G = 

Again we calculate the infinitesimals of the Lie symmetries (up to order 3) of the 
original system with coordinates {j/i, j/2, 2/3} and list them as X — {£, r)%, rj 2 , ^3}, 
representing the infinitesimal generator X = £dt + f]id yi + i] 2 d V2 + r}^d y3 . 

X 1 = {1,0,0,0} 

X 2 = {y 2 ,-k 1 y 1 y 2 ,k 1 y 1 y 2 ~£y 2 y 3 ,k 2 y 2 -£y 2 y 3 } 

X 3 = {2/3, -kmy3, hym - £y 2 yj,k 2 y 2 y 3 - £y 2 y 2 } 
k\ —ki —k 2 

x A = {0,— yi, — 2/1 + 2/22/3, — 2/2 + 2/22/3} 

£ £ £ 

,— 1 ki —k\ —k 2 
A 5 = { — t,—tyi, tyi+ty 2 y 3 , ty 2 + ty 2 y 3 \ 

£ £ £ £ 

— 1 k\ o —k\ 2 ~k 2 
X 6 = { 2/1,-2/1, 2/1+2/12/22/3, 2/12/2+2/12/22/3} 

£ £ £ £ 

Then we transform these to the u,v,w coordinate system as above with X = 
{£> rfi, ifc, V3} where X = £d t + ffid u + fj 2 d v + fj 5 d w : 

li - {1,0,0,0} 

X 2 = {v — w, — (v — w)n(u, v, w), — £(v — w)p(u, v, w), — ki(v — w)w} 
k\ 

~ f &2 1 . \ , \ i ^2 \ / \ 7 / ^2 \ 

A3 = {u + —w, — (u 4- — w)n(u, v,w) % —eiu + —w)p{u, v, w), — k\ (u + — it?) 

~ 1 

X4 = {0, — — n(u, «, ttfLpftt, v, w), — w\ 

£K\ £ 
—t —1 fci 

X5 = { — , — — tn(u, v, w),tp(u, v, w), — tw} 

£ £Kl £ 

Xg = { w, — — wn(u, v, w), wpiu, v, w), — w 2 } 

£ £Kl £ 

where n(u,v,w) — £k 2 w(—v + w) + k\{k 2 v + £u(—v + w)) and p(u,v,w) — 
(v —w)(u + j^w). 

5 The connections between the center manifold 
and the Lie symmetry 

Recently Cicogna and Gaeta [2] have written about the connections between 
dynamical systems and Lie symmetries. We are interested in particular in their 
results on invariant manifolds. They have commented that any Lie symmetry 
of the system will leave invariant both the stable and unstable manifolds. Due 
to the non-uniqueness of center manifolds, a Lie symmetry will map a center 
manifold to another (possibly the same) center manifold. The following result 
indicates when a center manifold given by uj(u) will be invariant under a given 
Lie symmetry, in their notation X = <j)d u + ipd v . 
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Lemma 5.1 (Lemma 4 of [5] Chapter 7). A center manifold w(u) is invariant 
under a Lie symmetry X = <f>d u + ipd v if and only if 

ip(u,u(u)) = (d u (uj(u))) ■ (f>(u,u)(u)). 

For the case G ^ 0, ui(u) = {0,0} giving zero on the right side of this 
equality. Thus the left side of this equation evaluated on the center manifold 
must always be zero if our center manifold is to be invariant under the action of 
the symmetry. This is the case with all of our Lie symmetries as given above. 
For example consider X 2 with <j>(u,v,w) = ^{G — k2)j{u,v,w)f(u,v 1 w) and 
ip(u, v, w) = {j(u, v, w){Gv + f(u, v, w)), —kiwj(u, v, w)}. Since j(u, 0, 0) = 0, 
tp(u,LLi(u)) — ip(u,0,0) = {0,0}, thus satisfying the necessary and sufficient 
condition of the lemma. It is easy to determine that the remainder of the 
symmetries in this case also leave the center manifold invariant. Thus the center 
manifolds inherit these Lie symmetries. However, in this case, since v = w = 0, 
all of our symmetries become trivial. 

Recall that in the case G = we found a center manifold to be the uv— plane. 
Now, in the notation of our lemma, lo(u) = 0, and again the right side of our 
equation is zero. Thus we must have ^(u^O) = for any symmetry that 
leaves invariant this center manifold. It can be easily checked to see that all 
of the symmetries listed above do indeed satisfy this necessary and sufficient 
condition. In this case the center manifold again inherits the Lie symmetries 
which are now nontrivial. The restriction of the system to our center manifold, 
the uv— plane, is 

v! = k2V — euv (27) 
v' = -euv. (28) 

The nontrivial symmetries inherited by this system are 

X 2 = {v, — v n(u, v, 0), — evp(u, v, 0), 0} 
X3 = {u, — un(u, v, 0), — eup(u, v, w), 0} 

h 

X 4 = {0,^n(u,«,0),i>(«,t;,Q),0} 
ek\ 

X 5 = {—t,^tn{u,v,0),tp(u,v,0),0} 

If we transform back to our original variables, we see that on the center 
manifold u = j/3 and v — y 2l resulting in the system: 



2/2 = -£2/2?/3 

2/3 = k 2V2 - £2/22/3 



(29) 
(30) 
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and the symmetries: 

X 2 = y 2 d t + {k 2 y\ - eyjy 3 ) d V2 - ey%y 3 d V3 (31) 
X 3 = yzdt + (k 2 y 2 y 3 - ey 2 yl) d V2 - ey 2 y\d V3 (32) 

x-a = (j-r-y? + viv^j d v2 + y2V3d y3 (33) 

X 5 = ~~ td t + (^T^ y2 + ^22/3^ d V2 + ty 2 y 3 d y3 (34) 

While we have calculated the infinitesimal generators, it would be interesting 
to determine the actual Lie symmetries on the center manifolds. We would like 
to say precisely what these maps do to various trajectories on the center manifold 
and to the flow in general. This is however, a very difficult question. There is 
no known method that allows us to take the infinitesimal generators of any Lie 
symmetry and integrate them to find the actual symmetries. The difficulty of 
this question is analogous to the solving of a system of differential equations 
analytically. 

For example, if we consider X 4 with r] 2 (t,y 2 ,y 3 ) = = ^y 2 + y 2 y 3 and 773 = 
2/22/3, this means that, letting 7 be the parameter of the one-parameter Lie 
group, we need to solve the following for y 2 and 2/3, giving us the map (2/2,2/3) 
as our symmetry: 

dy 2 -k 2 , _ . , 

-3— = 2/2 + 2/22/3 (35) 

d'y e 

dy 3 . „ , , 

d7 = V2 V z (36) 

This is equivalent to the system above. Attempting to solve this system we find 
it equivalent to solving the following: 

2/2 = e^-^+y^ (37) 
2/3 = ef*"* (38) 

with the initial conditions 2/2(7, y 2 , y 3 )\ 7 =o = 2/2 and 2/3(7, y 2 , 2/ 3 )| 7 =o = 2/3- 

This is something we continue to work on for this particular system as well 
as in general. 



5.1 Comparison to previous results 

Based on numerical solutions of the original system of equations Ross et. al 
[3 predicted trajectories for M, M*,A and D with particular emphasis on the 
concentrations of M* (cells undergoing multiplication) and A (the antagonist). 
They found that the behavior depended on the values of the various constants 
fcj. In particular, with k 3 = and G > 0, they found unrestrained growth of 
both M* and A. For the values k 3 = Q and a negative G, A increases toward an 
upper limit and M* increases slightly but then begins to decrease toward zero. 
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For &3 > and G > 0, both M* and A increase initially, but then M* reaches 
a maximum and begins to decline while A approaches an upper bound. All of 
these analyses combined to indicate to the food scientists that the necessary 
constraints for growth-death kinetics are non-zero values for ks and positive 
values of G. 

In our consideration of the system, we also found that M*(= y^) and A(= y^) 
were the two variables that determined the behavior of the system. In the G^O 
case, the center manifold is the u-axis, which corresponds to A when all other 
variables are zero, as on the center manifold. When G = 0, the reduced system 
on the center manifold is given by equations [27] and [28j An inspection of this 
system, noting that u = A and v = M* , shows that the behavior is qualitatively 
identical to that found numerically in [4] for the case k = [1 4 100 4], i.e. 
k-2 = ki — 4 resulting in G = 0. In both the results are that M* goes to zero 
and A approaches a constant value. 
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